AN ANALYSIS OF NODE-BASED CLUSTER SUMMATION 
RULES IN THE QUASICONTINUUM METHOD 



MITCHELL LUSKIN AND CHRISTOPH ORTNER 

Abstract. We investigate two examples of node-based cluster summa- 
tion rules that have been proposed for the quasicontinuum method: a 
force-based approach (Knap & Ortiz, J. Mech. Phys. Solids 49, 2001), 
and an energy-based approach which is a generalization of the non-local 
quasicontinuum method (Eidel & Stukowski, J. Mech. Phys. Solids, to 
appear). We show that, even for the case of nearest neighbour interac- 
tion in a one-dimensional periodic chain, both of these approaches cre- 
ate large errors when used with graded and, more generally, non-smooth 
meshes. These errors cannot be removed by increasing the cluster size. 
We offer some suggestions how the accuracy of (cluster) summation rules 
may be improved. 



1. Introduction 

The quasicontinuum (QC) method [15, 16, 18, 19] is a prototypical coarse- 
graining technique for the static and quasi-static simulation of crystalline 
solids. One of its key features is that, instead of coupling an atomistic model 
to a continuum model, it uses the atomistic model also in the continuum 
region where degrees of freedom are removed from the model by means of 
piecewise linear interpolation. 

However, the nonlocal nature of the atomistic interactions makes further 
approximation necessary to enable the computation of energies or forces 
with complexity proportional to the number of coarse degrees of freedom. 
Two families of approximations have been developed to achieve this goal. 
One family of approximations localizes the interactions by a strain energy 
density (based on the Cauchy-Born rule) which provides sufficient accuracy 
in regions away from defects where the strain gradient varies slowly. Classi- 
cal finite element methodology can then be utilized in those regions modeled 
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by the strain energy density. This class of quasicontinuum approximations 
have been the subject of many recent mathematical analyses [2-7,12,13,17]. 

The purpose of the present paper is to investigate the second family of 
approximations that have been developed to reduce the computational com- 
plexity of the quasicontinuum method. These methods, which have received 
far less mathematical attention, use summation rules (discrete variants of 
continuous quadrature rules) to approximate the sums that define the qua- 
sicontinuum energy or forces. To the best of our knowledge, the force-based 
cluster summation rule of Knap and Ortiz [11], the non-local QC method 
based on a simple trapezoidal rule [15, Sec. 3.3], or its extension to energy- 
based cluster summation rules [8] have not been analyzed to date. These 
cluster summation rules approximate the sum over atom-based quantities 
by uniformly averaging over atoms in clusters (balls) around the nodes and 
then by weighting these cluster averages so that summands that are obtained 
from piecewise linear interpolation with respect to the quasicontinuum mesh 
are exactly computed. 

In a recent benchmark of different QC methods [14], the cluster summa- 
tion rules do not compare favorably with quasicontinuum approximations 
that utilize the strain energy density or with other atomistic-to-continuum 
coupling methods. In the present paper, we give a simple yet rigorous ana- 
lytical explanation for this poor performance. We demonstrate that, even for 
the simplest imaginable atomistic model, a periodic one-dimensional chain 
with harmonic nearest-neighbour interaction, the cluster summation rules 
formulated in [8,11] lead to inconsistent and inaccurate QC methods when 
used with graded and non-smooth meshes. Increasing the cluster size does 
not resolve this problem. The benchmark [14] uses a mesh that is refined 
from large triangles to the atomistic scale as is typical for quasicontinuum 
computations. Our analysis shows that this kind of refined mesh would 
lead to an inconsistent and inaccurate method for the approximation of our 
one-dimensional model by cluster summation rules. 

The atomistic model, the finite element space {coarse space), and some 
additional notation are introduced in Sections 1 1.11 and 1 1.21 We treat the two 
classes of cluster summation methods separately; the force-based summation 
rule in Section [2] and the energy-based summation rule in Section [3l These 
sections can be read independently of each other. Finally, in the conclusion, 
several possibilities are identified how cluster summation rules might be 
modified in order to lead to accuracte QC methods. 

Although our analysis treats the approximation of the discrete sums in 
the quasicontinuum energy, the reasons for the lack of accuracy of the clus- 
ter summation rules can already be understood from applying the cluster- 
method concepts to the finite element approximation of continuum elastic- 
ity [1]. The force conjugate to a finite element nodal degree of freedom (the 
negative of the partial derivative of the elastic energy functional constrained 
to the space of finite element trial functions) depends only on integrals of the 
jump of the displacement gradient across the element boundaries (the finite 
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element force density). The accuracy of cluster-based quadrature relies on 
the finite element force density being smooth in space; however, the support 
of this force is concentrated on the element boundaries. Thus, the conjugate 
forces obtained from the cluster-based quadrature will be much too large. 

Accurate node-based quadrature rules used in the approximation of a 
continuum finite element energy are computed within each element during 
the assembly process, which leads naturally to an accurate mesh-dependent 
non-uniform weighting of the energy density in any ball surrounding each 
node. The cluster-based quadrature approximation uses a uniform weighting 
in the ball surrounding each node. Thus, since the energy density for a 
displacement in the finite element trial space is generally discontinuous at 
the nodes, the cluster-based quadrature rules are likely to be inaccurate for 
nonuniform meshes. 



Remark 1. We have not included the formulations of Lin [13] or of 
Gunzburger & Zhang [9, 10], which are closely related to cluster summa- 
tion methods, in our analysis. Our main reason for this exclusion is that 
their formulations do not suffer from the same deficiencies as the methods 
which we investiate here. If errors are present in their approach (we have 
not investigated this further), they would most likely be caused by finite 
range interaction and cannot be observed for the simple nearest-neighbour 
interaction system which we investigate here. 

1.1. The model problem. We choose the simplest immaginable atomistic 
model problem, a one-dimensional periodic chain with nearest-neighbour 
pair potential interaction. The continuum reference domain will be (—1,1]. 
For fixed N € N, the atomic spacing is given by e = 1/A^, and the atomistic 
reference lattice by 



We refer to Remark [2] for a motivation of the constraint vq = 0. 

We assume that the interatomic interaction reaches only nearest neigh- 
bours, and that the only external force is a dead load. Thus, we can write 
the total energy as a sum of a stored energy £{v) and an external potential 
energy -f[v\ : 



C = {e£:e = -N + 1,...,N}. 



The space of periodic displacements of C is denoted 



X = {v : ve+2N = vi for £ £ Z, and vq = O}. 





(2) 
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where 



£{v) 



N 

e(l){e'^ {vi - Vi-i)) , and 

e=-N+i 



N 



m 



-N+l 



Here (/> is assumed to be smooth, at least in a neighbourhood of zero, and 
{fi)eez is a fixed 2A^-periodic sequence. In fact, we shall usually make the 
simplifying assumption that (p is a convex quadratic and that / is obtained, 
for example, by nodal interpolation from a smooth 2-periodic function f{x). 
Such assumptions are valid for small displacements from the reference state. 

Rescaling the domain (and the energy) by the atomistic spacing e is not 
strictly necessary, but it helps us understand the connection of the atomistic 
problem to continuum theory. 

The atomistic problem is to find 



where 'argmin denotes the set of local minimizers of $ in X. The first 

order necessary criticality condition, in variational form, is 



in short, £'{u) = f. 

Remark 2. In the definition of the displacement space ([T]), we have 
imposed the condition vq = for admissible displacements. This is one of 
several ways to remove the zero mode from the space in order to render 
the energy functional <I> coercive. Furthermore, this constraint allows us 
to easily construct a problem with a 'singularity' at the origin (cf. Section 
13. 2p . In general, if the external force / is 'smooth' and anti-symmetric, then 
the solution will be 'smooth' as well. If / is not anti-symmetric, then the 
solution may have a 'kink' at the origin even if / is smooth. 

We fix some additional notation, some of which we have already used 
above. The arguments of nonlinear functionals are enclosed in round brack- 
ets while those of (multi-)linear forms are enclosed in square brackets, for 
example, £{u) or f[u]. The Frechet derivatives are denoted by ', for exam- 
ple, S'{u) is a linear form on X. Consequently, £'{u)[v] denotes a direc- 
tional derivative. Similarly, £"{u) is a bilinear form on X , and it is written 
£"{u)[v,vj\ with arguments v,w S X. Finally, we will frequently use the 
notation v[ = e~^{v£ — Vi^i) to denote the differences. 

Atomistic displacements are always identified with their piecewise affine 
interpolants. In particular, for v G X,we have v'^ = v'{x) for x G {{( — 1)^, ^e) 
Through this identification, the space X is naturally embedded in the spaces 



W\f{-1, l) = {vG W^'P(M) : v{0) = 0, v{x + 2) = v{x) Vx E R}, 



u G argmin ^{X) 



(3) 



£'{u)[v]= f[v] yv£X, 



(4) 
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for 1 < p < oo. 

1.2. The constrained approximation. The quasicontinuum approxima- 
tion to the atomistic model problem ([3]) is obtained in two steps: (i) replacing 
the displacement space X hy a low-dimensional coarse space X^, and (ii) ap- 
proximating the nonlinear system for arguments from the coarse space. 
Often, the process is in fact reversed, however, for the class of QC methods 
which we consider in the present paper, the order is as stated above. 
We fix a set of rep-atoms 

Ch = {eik--k = -K + l,...,K} cC, 

so that #Ch <C #£. The set is 2K-periodically extended, that is, we 
define lk+2K = + 2A^ for all A; € Z. The coarse space can therefore be 
written as 

'^/i = {'W/i £ : ^^/i is piecewise affine with respect to (e^fc)fcGz}- 

The constrained atomistic approximation is to find 

Uh e argmin$(A'/i), (5) 

for which the first order criticality condition is 

£'{uh)[vh] = f[vh] G Xh. (6) 

Even though the number of degrees of freedom is significantly reduced in ([5]) , 
the nonlinear system ^ is still prohibitively expensive to evaluate since it re- 
quires summation over all atoms. Hence, the second step of the QC method, 
the approximation of the nonlinear system, is as important as the coarsening 
step. One class of methods to achieve this are the (cluster-)summation rules 
which we investigate in the following sections [8,11]. 

We conclude the introduction with some additional notation related to the 
coarse space X^- For k £ Z, we denote = e{£k — ^fc-i)- We will assume 
throughout that the mesh satisfies a local regularity condition: there exists 
a constant k > 1 such that 

i^^^hk-i ^ hk ^ i^hk-i for k = —K -|- 1, . . . , (7) 

For Vh e Xh, we denote V = {Vk)kei = {vejkez the (2i^-periodic) vector of 
nodal values, so that 

K 

Vh,e= Yl ^kCkiei), i=-N + l,...,N, (8) 

k=-K+l 

where Cfc denotes the periodic nodal basis function associated with node 
e£k- Furthermore, we denote = (14 — Vk-i)/hk the gradient of Vh in the 
element (e£fc-i, e^fc)- In particular, we have v'f^g = Vj^ if £k-i < i < ik- 
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2. FORCE-BASED SUMMATION RULES 

Since they are more easily understood, we shall first investigate the force- 
based summation rules, introduced by Knap and Ortiz [11]. Our presenta- 
tion closely follows their formulation. 

Instead of viewing the constrained approximation ([5]) as a minimization 
problem, we concentrate purely on the equilibrium equations ([6]). However, 
instead of the variational formulation ^ , we use the nodal force formulation 

Using the expansion ([8]) of Uh in the nodal basis (Cfc)fcL-ic'+i' nodal 
forces are rewritten in the form 



, u=uu dU i 



dUj „ ^ _ duf 



dU-i , ^ du 



that is, 



^= E F,iu,)CM). (10) 







where 



FAu)- ^ ' 



dui 

At this point, we apply a cluster summation rule to approximate the sum 
in (dOD, 

K 

Fj,h{uh)-= J2 ''kY.F.iuhKjiei), j = -K + l,...,K, (11) 

where the sets Ck are clusters surrounding the repatoms £k, 

Cfc = {4 - r^, . . . ,4 + r+}, k = -K + l,...,K, 

and the weights z/^ are defined by the requirement that the basis functions 
are summed exactly. 



N K 

l=-N+l k=-K+l 



Cj{e£)= J2 j = -K + l,...,K. (12) 



In practice, the system (jl2p is solving using a mass lumping approximation 
[11, Sec. 3.2] yielding approximate weights z^^. We will only investigate the 
effect of the cluster summation rule in the situation when e <C /ifc for all k, 
hence we shall assume throughout that = r for all k. In this particular 
case, we show in Appendix I A. II that 

" 2(2r + 1)£ ' = ^k + 0(1). 
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We note that the relative error is of order 0{re/[hk + /ifc+i))- 

2.1. Analysis without external forces. In this section, we assume that 
/ = 0. To motivate this assumption, we note that external body forces are 
only rarely the driving force in an atomistic simulation and would therefore 
distort the picture we are about to present. We shall simply ignore the fact 
that, as a result, the atomistic problem becomes trivial. 
If / = 0, it follows that 

F,{u) = <P'{u',) - <P'{u',^,), e = -N + l,...,N, 

and, if we insert u = £ X^, we obtain 

It follows that, independently of the cluster size, we obtain 

FkM^h) = i.k{<P'{Ul,)-cP'{U;,^,)). 
Since ^ 0, for all k, the equation Fk^h{uh) = is equivalent to 

Qjj^ = <t> [Uk) - <t> (f^fe+i) = 0. 

Thus, we see that, even though the cluster summation rule (jlip is grossly 
inaccurate for moderate cluster radii (the weights are of order 0{hk/{re))), 
the resulting system is nevertheless equivalent to an exact evaluation of the 
full constrained approximation which is known to be an excellent approxi- 
mation to the full atomistic system [17]. 

Remark 3. We need to be careful in extrapolating this observation to 
the case of finite range interaction and indeed the much more subtle and 
interesting two- and three-dimensional setting. These situations need to be 
investigated in more detail. Nevertheless, we can make some comments to 
motivate further investigation. The main observation that we have made 
in the present section, that forces are concentrated on the interfaces is still 
valid. In 2D and 3D, and for general atomistic models, the identity 



dUk ^ due 

remains true (note, however, that now C C M'^, d € {2,3}, and u : C ^ 
Furthermore, in the absence of an external force, in the interior of 



a large element, the force ^q^^ \ u=Uh is zero (or negligibly small) for an 
atomistic model with short range interaction. From this, we see that the 
contributions to the force on the node ik are concentrated near all element 
faces which touch the repatom £k ■ This shows that the summation rule used 
to obtain Fk^h{uh) should be obtained from a summation over these faces 
(surface integration) instead of summation over the entire patch (volume 
integration) . 
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2.2. Analysis with non-zero forces. We now return to the case of non- 
zero forces. We shall assume that the forces {fe)t(^z are obtained by in- 
terpolating a smooth 2-periodic function / G C^[— 1,1]. In this case, we 
obtain 

F,{u) = - -efi, £ = -A^ + 1, . . . , AT, 

and hence, 

F,K) = |^'(^^)-^'(^^-^)-^^' 'T''^ (14) 
^ ^ \ -efi> , otherwise. ^ ' 

It is then fairly straightforward to see that 

Fk,h{uh) = T^k{4>'{K) - (P'iK+i)) - fk, 

where fk is obtained by applying the cluster summation rule to the external 
forces only, 

K 

fj= E '^kY.eMjiel), j = -K+l,...,K. 
k=-K+i ^eCfe 

Since we assumed that / G C^[— 1, 1], we can deduce from a fairly straight- 
forward interpolation error analysis (cf. Appendix IA.2P that 

N 

h= eM,{ei) + 0{h] + h]^^) = f[C,] + 0{h] + h]+^). (15) 

l=-N+l 

The force-based cluster quasicontinuum equations (llip therefore become 
uu{<t^'{U'^) - <A'(f/fc+i)) = /[Cd + 0{hl + hl^,), 

as opposed to the 'exact' equations of the constrained quasicontinuum ap- 
proximation 

To illustrate this point further, let us assume that the interaction is har- 
monic, that is (j){r) = ^r^, and that the mesh is uniform (/i^ = h for all 
k). In that case, the weights are given by Uk = h/{e{2r + 1)) (cf. Appendix 
IA.1|) . and hence 

Since the difference operator on the left-hand side is linear, we therefore 
deduce that 

n 

where u/j is the solution of the constrained atomistic system ([6]). In the 
typical case when er <C /i, this result demonstrates the catastrophic error 
made in the force-based cluster summation rule. The reason why we do 
not observe a similar cancellation effect as in Section 12.11 is because the 
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external force contribution was summed accurately, while the summation of 
the interatomic forces is grossly inaccurate. 

3. Energy-based summation rules 

We have seen in the previous section that the failure of the cluster sum- 
mation rule applied to the force balance equations fails because a 'volume 
integration' method was applied to a 'surface integral'. It is natural, there- 
fore, to investigate the cluster summation rule applied to the energy func- 
tional. This would lead to a conservative coarse grained system, which was 
the main motivation for Eidel and Stukowski [8] to use this method. They 
have noted in [8, Sec. 5] that this method also has shortcomings, and we 
shall analyze these in detail in the present section. 

To formulate the energy-based cluster summation rule, we first rewrite 
the stored energy functional E in the form 

N 

£=-Af+l 

where 

The term Eiiu) is the contribution of the atom at site i to the overall energy. 
The sum over the terms Et{u) is approximated by a summation rule of the 
form 

N K 
£=-Ar+l k=-K^\ l&C^ 

where the sets = {^^ — r^, . . . + r^} are non-overlapping clusters 
surrounding the repatoms. The weights uJk are determined by requiring 
that the summation rule is exact for all basis functions, that is, 

N K 

J2 eC,{e£)= ^kY^Cjiei), j = -K + 1, . . . , K. (17) 

£=-N+l k=~K+l l&Ck 

To motivate a simplification which we are about to make, assume, for the 
El] that 

hk + hk+i 



moment, that = r for all k. For this case, we have shown in Appendix 



Furthermore, we observe that 

^ £t{vh) = rct>{Vi) + \{ck{yi) + </.(yfeVi)) + rct^iVUi) 

= i(2r + l)(</.(yfc')+<^(^m)) 



= {2r + l)8,^{vh). 
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Thus, we see that, up to an error of order 0{e), a finite cluster size reduces 
immediately to a discrete trapezoidal rule. 

In view of this observation, we shall assume throughout this section that 
Ck = {^k}- The approximate energy functional becomes 

K 

£h{vh) = ^ UkSlki'^h), (19) 
k=-K+l 

with weights uok = ^{hk + hk+i)- This method (fT9]) is sometimes labelled 
the non-local QC method [15, Sec. 3.3]. 

Remark 4. 1. The observations made above are only partially valid for 
non-nearest neighbour interaction. In that case, additional interface terms 
of the form 4>{V^ + V^_^^) enter the QC energy functional. 

2. A further correction from our simplifying assumption needs to be taken 
into account when the mesh is refined to atomistic level where we need to 
use variable cluster sizes. For simplicity, we have chosen to ignore this 
further complication, but note that our analysis in Appendix lA.ll can be 
generalized to variable cluster sizes provided the cluster radii are symmetric 
in each element (that is, r^_i = t^). In that case, we would obtain a similar 
formula as (llSp but with an error of order O(emaxfcr^). 

For simplicity, we assume that the dead load f[v] is not approximated. 
The total energy for the QC method is therefore given by 

^h{vh) = Shivh) - f[vh], 
where is defined in (jl9p . The corresponding criticality condition is 

£hiuh)[vh] = f[vh] yvh£Xh. (20) 

In order to analyze the (non-local) QC method, we first rewrite £h in the 
form 

K K 

£h{vh)= ^kW{Vi) + <l^{VU))= \{^k+^k-MVl). 

k=-K+l k=-K+l 

Using uJk = ^{hk + hk+i), we see that 

liujk + Wfc-i) = lihk-i + hk+ h + hk+i) 
= hk + \{hk^i - 2hk + /ife+i) 
= : hk{l + u)fc). 



where 

and hence, we obtain 



hk-i — + hk+i 

^k = 77 , (2i) 



K 



£h{vh) = £{vh) + hk^kW). (22) 



fc=i 
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We use uj to denote the piecewise constant function taking values ilik in the 
elements {e£k-i, £^k)- 

We can relate the connection between the error in the cluster approxi- 
mation to the error in the trapezoidal rule by noting that (|22p is equivalent 
to 

K 

£h{vh) = £{vh) + i {<P{vUi) - mvi) + <t>{yLi)) ■ (23) 

k=l 

Prom (|22|) and (j23|) , we already anticipate that the local mesh smoothness 
will have a significant impact on the accuracy of the method. For example, 
if the cifc are oscillatory, then it is possible to lower the energy by introducing 
a 'microstructure' into the quasi-continuum displacement. We will see this 
in detail in Example 3 below. In Example 2, we discuss another effect that 
may introduce large errors in the simulation. 



Remark 5. We could analyze the consistency of the method using finite 
difference techniques. Taking the derivative of £h{uh) with respect to the 
nodal value f/^, we obtain 

where lu^ is given by (j21|) . One can see here as well that, if the terms cDfc 
are not close to zero, then the method is inconsistent. However, a rigorous 
error analysis is more conveniently performed in the variational setting of 
the finite element method. 

To further simplify the analysis, we assume from now on that the inter- 
action is harmonic, that is, (f){r) = ^r"^. This assumption can be justified, 
for example, for small perturbations from a reference state. In that case, 
the fully atomistic problem ([4]) has a unique solution u. Furthermore, let 
Uh be the unique solution of the constrained approximation, which is the 
best approximation to u from in the energy norm. Since all weights w^, 
k = 1, . . . ,K, are positive, it follows that the QC functional <I>/j also has a 
unique critical point, Uh- 

Since £ and £h are both quadratic, their Hessians are independent of the 
arguments. Thus, we will write £"f^s^[vh,Wh] instead of, say, £"{uh)[vh,'Wh]- 
With this notation, the criticality conditions ([6]) and ()20p become, respec- 
tively, 

£"[uh,vh] = f[vh\ e Xh, 

£h[uh,Vh] = f[vh] Vu/j E Xh. 

Thus, the error Uh — Uh satisfies, for all Vh G X^, 

£h[uh - Uh, Vh] = f[vh] - S'hiuh, Vh] = £"[uh, Vh] - £h[uh,Vh]. 
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Using Strang's First Lemma [1, Thm. 4.1.1] and the mesh regularity condi- 
tion ([7]), we obtain 

5(1 + - u'^llLa < {S'jlluh - Uh,Uh - Uh])^^"^ (24) 

< sup \£"[uh,Vh] - £h[uh,Vh]\. 

We wish to obtain sharp upper and lower bounds on the variational crime. 
To this end, we first note that a piecewise constant function r/j, is the gra- 
dient of an element of Xh if, and only if, J^-^^rhdx = 0. Thus, setting 



K _ /-i 



(25) 



k=-K+l 

we obtain 

sup \£"{uh,Vh) - £h{uh,Vh)\= sup / [ujuy^- a)v't^Ax 

= lli^^h - a||L2 =: p{uh), (26) 

which gives an upper bound on the error. To obtain a lower bound as well, 
we reverse the argument in ()24p . yielding 

p{uh) = sup \£"[uh,Vh]- £h[uh,Vh]\ 

= sup \£'/^[uh-Uh,Vh\\ 

< ^{1 + ^^)\\u'h-u'h\\L2■ 
Comb'mmg this estimate with (|24p and (j26p . we finally arrive at 

i(l + K-i)||4 - <||l2 < Kn/,) < 1(1 + k)||4 - <||l2, (27) 

where 

K 

fc=-_ft'-i-i 

Note that (j27p does not estimate the actual error u — but the deviation 
from the best approximation in the energy norm. In the following, we will 
investigate the term p{uh) for three typical meshes that may occur in a 
simulation. 

3.1. Example 1: smooth meshes. We begin by looking at a somewhat 
idealistic situation. We assume that e <^ h and that the mesh nodes at ei^, 
k = —K -|- 1, . . . , /C, are given by a smooth (and periodic) deformation </? of 
the periodic domain (—1, 1], that is, 

e£k = (f{hk), k = -K + 1,...,K, 
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where h = 1/K, and fix) = (p{x) + 2 for all x G M. In that case, the term 
Uk can be estimated, using Taylor's Theorem, to obtain 

hk+i — ^hk + hk-i 



<p{h{k + 1)) - 3(p{hk) + 3(p{h{k - 1)) - ip{h{k - 2)) 



Aiipihk) - ifiih - l)k)) 

where Xk = {k — ^)h. In particular, we obtain 
where C = |maxj.g[_i ij \(p"'{x)/ip'{x)\. Since 

K 

P{uhf = XI ^k\'^kU'k\ -2a^, 
k=-K+\ 

it follows that 

K 

p{uhf< X hk\(bkU'k\^ <c''h''\K\\l,+o{h!'), 

k=-K+l 

which is of a smaller order than the best approximation error. 

3.2. Example 2: graded meshes. Since the main target of the quasicon- 
tinuum method are problems with defects or singularities, extremely smooth 
meshes satisfying the assumptions of Example 1 are rare in quasicontinuum 
applications. The most important example for the QC method is a mesh 
which refines to atomistic level. To investigate this situation, we construct 
an exponentially graded mesh as follows. We fix iiT > and N = 2^~^, and 
define = and 

4 = sgn(A;)2l'=l-\ k = -K + 1, . . . , K. 
In that case, we obtain 

£, k = 0, 1, 
hk=\ 2'=-2£, k = 2,...,K, 

_ 2l*^l-ie, A: = -l,-2,...,-i<: + l, 

which, in particular, gives the mesh regularity parameter n = 2. We can, 
furthermore, explicitly compute the coefficients 

0, k = 0,l, 

1/4, k = -l,2, 

1/8, A: = -i^ + 2,...,-2,3,...,i^-l, 

-1/8, k = -K+l,K. 

To further investigate the error, let us assume that the displacement gra- 
dient in the 'continuum region' is negligable. Let us further assume that 



i^k 
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-1 -0.5 0.5 1 

Figure 1 . Computational example on a highly graded mesh 
with force f{x) = 10^ exp(-10''x2), N = 2^4, and K = 15. 
Since / is not anti-symmetric, the solution is non-smooth 
at the origin (cf. Remark [2]) . The relative error in the en- 
ergy norm satisfies — "u'^Hl^/II^/iIIl^ ~ 0.11 (that is 11%), 
which is in excellent agreement with our prediction. The rel- 
ative error for the energy satisfies {£{u) — £h{uh))/\S{u)\ ^ 
—0.13 (that is 13%). Precisely as we predicted, we see that 
the failure in the coefficients enforces a smaller QC displace- 
ment which results in a higher energy. 

the displacement gradient does not vary considerably between the elements 
(0, h) and {h, 2h) as well as between {—h, 0) and {—2h, —h). In that case, we 
can ignore the Cok = >^-K+i = — | coefficients in the outmost elements, and 
we can 'split' the coefficients uj2 = among the purely atomistic elements 
in order to obtain a ~ and 

K K 

p{uhf = ^ hk\(l:kU'k - a] hk\tJ'k\ = {\\\uh\\i,2) . 

k=-K+l -K+l 

Prom (I27p . we would therefore expect a relative error of size 1/8 (that is 
12%), independent of the mesh size h. This is in perfect agreement with the 
computational example we present in Pigure [H 

3.3. Example 3: a non-smooth mesh. In the final example of our anal- 
ysis of the energy-based cluster summation rule, we consider a mesh which 
is quasi-uniform, but not smooth. We assume that e <^ h, and that 

/ifc = i/i(3 + (-l)'=), k = -K + 1,...,K, (28) 




ANALYSIS OF NODE-BASED CLUSTER METHODS IN THE QC METHOD 15 



that is, we have hi = h, h2 = 2h, = h, and so forth. A mesh of precisely 
this type wih rarely be found in practise, however, it is an excellent model 
situation that demonstrates a source of error for non-smooth meshes. 
In this situation, the coeffcients tOk satisfy 

-1/4, if A; is even, 
1/2, if A; is odd. 

Suppose that Uh is the interpolant of a smooth function u, so that U'f, varies 
little between elements. Then the oscillatory nature of the coefficients Uk, 
weighted according to the size of the elements, indicates that a ~ 0. More 
precisely, we have 

i(/ifcC7^ + /ifc+iC/^+i) + lhk+i{UUi - U'k). k odd, 
-|(/^A:-iC/fc_i + /ifc^fc) + IhkiU'j, - t7^_i), k even, 

from which we can easily deduce that \a\ < Ch, where C depends on the 
second differences of (which we assumed to be moderately small). Some 
similar, algebraic manipulations show that 
K 



hk'^kUl 



E 



k=~K+l 

and thus, we obtain 

K 

k=-K+l 

From our general error estimate (j27p . we therefore expect the relative 
error to be roughly of the order l/\/8 (that is 35%), which is in excellent 
agreement with the numerical example shown in Figure [2l 

4. Conclusion 

We have shown that node-based cluster summation rules, applied either 
to the force-based formulation of the QC method or to the energy-based for- 
mulation of the QC method lead to inconsistent and inaccurate numerical 
schemes when used with graded or non-smooth meshes. We stress, further- 
more, that increasing the cluster size is not a remedy for the sources of error 
which we have discussed. 

We do not rule out, however, that QC methods based on a more careful 
choice of summation points may yet lead to excellent computational tools. 
We would like to comment on three options which qualify for further inves- 
tigation: 

Lin's formulation [13] and the formulation of Gunzburger &; Zhang [9,10], 
which are based on summation points in the interior of the elements do 
not suffer from any of the deficiencies which we have found in the present 
work. It will be necessary, however, to carefully investigate the effect of 
next-nearest neighbour and finite range interaction in the transition region 
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Example 3 




Figure 2. Computational example on an oscillatory mesh 
with force f{x) = sin(7rx), N = 10"^ and K = 20. The fully 
atomistic solution is given by u{x) = tt"^ sin(7ra:). The rela- 
tive error in the energy norm satisfies — 'S^||L2/||^i/j||L2 ~ 
0.33 (that is 33%) while the relative error for the energy sat- 
isfies {£{u) - £h{uh))/\£{u)\ « 0.097 (that is 9.7%). In the 
zoomed box, we see that the microstructure, induced by the 
oscillatory coefficients, is lowering the energy and creates a 
'non-smooth' quasicontinuum solution. 



in which the mesh is refined from large triangles to the atomistic scale where 
all degrees of freedom are retained. 

A force-based formulation, where the summation is performed over el- 
ement interfaces rather than elements may yet lead to an accurate QC 
method. This is clearly true in one dimension, but needs to be carefully 
studied in higher dimensions. 

Finally, we propose to investigate the possibility of assigning variable 
weights to atoms within the same cluster, in an energy-based cluster sum- 
mation rule. It can be readily verified, following the analogy with continuum 
finite element energies discussed at the end of the introduction, that if the av- 
erage over atoms in a cluster is weighted according to element sizes, then the 
resulting method will be accurate for nearest-neighbour interaction. Once 
again, the crucial questions are whether this accuracy can be retained for 
finite range interaction and the application relevant two- three-dimensional 
situations. 



Appendix A. Proofs 



ANALYSIS OF NODE-BASED CLUSTER METHODS IN THE QC METHOD 17 

A.l. Computation of summation weights. The analysis presented in 
this appendix applies to both the energy-based and the force-based summa- 
tion rules, since the weights satisfy Uj = evj. For no particular reason, we 
chose to work with the weights {'^j)f=_K+v 

We assume throughout that = r. According to the requirement (jl7l) . 
the governing equations for the weights uj = {uJk)^=-K+i ~ 9 where 

K 

k=-K+i eeCk 
= Wj-i(^r(r + l)ehj^) + u;j+i(|r(r + l)ehjl-^) 

+ cjj((2r + 1) - ir(r + l)ehj^ - ir(r + l)ehjl^) 

and 

N 

e=-N+i 

To prove that M is invertible, we show that it is row-diagonally dominant. 
For each j, we have 

Mjj - \Mj,k\ = (2r + 1) - r(r + l)eihj' + hjl,). 

Since we assumed that the clusters do not overlap, it follows that £j — £j-i > 
2r + 1, in particular, —er > —\hj, from which we deduce that 

- Y \Mj,k\ > (2r + 1) - (r + 1) = r. (29) 

Thus, M is invertible and uj is well-defined. We note, furthermore, that (I29p 
implies that 

iiM-'iu»s{ Si: im 

Our next observation is that the lumped system for computing the ap- 
proximate weights (D = {ij^k)f=^K+i 

(2r + = + hj), j = -K + l,..., K, 

that is. 

We shall now prove that the exact weights are only 0{e) 

perturbations from the approximate weights obtained by mass-lumping. To 
this end, we define the residual p = {pk)k=-K+i' 

Pj ■= {Mu})j -gj 
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If the mesh is uniform or if r = 0, then p = 0. In general, under the mesh 
regularity assumption ^ we obtain the residual estimate 

< C{K)er. (32) 

To estimate the error on the weights, we note that M(ui — uj) = p, and hence 

||w — wll^oo < ||M~"^ ll^oo IIpII^oo < ma.x{l,r)^^C{K)re, 
that is, we obtain 

II- II ^ / C{K)e if r > 1, . . 

Il'^-'^ll.-<| ^0 ifr = 0. ^^^^ 

This may seem an impossibly strong result at first glance, however, we note 
that it is only true under the restriction that ^ r for all k. We expect 
that, in general, the cluster size will influence the estimate to give an error 
of order ©(emaxr^). 

Upon noticing that the weights for the force-based summation rule satisfy 
Uj = ujj/e, an estimate for the weights vj follows immediately. 



A. 2. Proof of (lisp . Let I be the interpolation operator for the quasicon- 
tinuum mesh, that is, I : X ^ with [Ivh)tf. = Vh/k^ ^ ~ ~^ + I, . . . , K. 
Let g £ X he given by g£ = ffC,j{et}, then 

N N N 

e=-N+l -N+l -N+l 

In view of our assumption that = f{e£), where / G C^[— 1,1], and the 
interpolation error estimate of [17, Thm A. 4], it follows that 

N N 

e=^N+i -N+i 

Since, by definition, piecewise affine functions are summed exactly by the 
cluster summation rule, it follows that 

N K 

egi= VkY,^l9i + 0{h] + h]^{). 

l=-N+l k=-K+l l€Ck 

Applying the same argument as above, we can deduce that 

K K 

k=-K+i eeCk k=-K+i eeCk 

from which the desired result follows. 
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